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Abstract 

We propose to extend the well-known MUSCL-Hancock scheme for Euler equa- 
tions to the induction equation modeling the magnetic field evolution in kinematic 
dynamo problems. The scheme is based on an integral form of the underlying con- 
servation law which, in our formulation, results in a "finite-surface" scheme for the 
induction equation. This naturally leads to the well-known "constrained transport" 
method, with additional continuity requirement on the magnetic field representa- 
tion. The second ingredient in the MUSCL scheme is the predictor step that ensures 
second order accuracy both in space and time. We explore specific constraints that 
the mathematical properties of the induction equations place on this predictor step, 
showing that three possible variants can be considered. We show that the most ag- 
gressive formulations (referred to as C-MUSCL and U- MUSCL) reach the same level 
of accuracy as the other one (referred to as Runge-Kutta) , at a lower computational 
cost. More interestingly, these two schemes are compatible with the Adaptive Mesh 
Refinement (AMR) framework. It has been implemented in the AMR code RAM- 
SES. It offers a novel and efficient implementation of a second order scheme for the 
induction equation. We have tested it by solving two kinematic dynamo problems in 
the low diffusion limit. The construction of this scheme for the induction equation 
constitutes a step towards solving the full MHD set of equations using an extension 
of our current methodology. 

Key words: 76W05 Magnetohydro dynamics and electrohydro dynamics, 85A30 
Hydrodynamic and hydromagnetic problems, 65M06 Finite difference methods. 
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1 Introduction 



The extension of Godunov-type conservative schemes for Euler equations of 
fluid dynamics (Toro, 1999; Bouchut, 2005) to the system of ideal magneto- 
hydrodynamics (MHD) has been a matter of intensive research, starting from 
the early 90's. The great variety of different MHD implementations of the 
original Godunov method, especially in a multidimensional setting, has left 
several unexplored paths opened in designing MHD conservative methods. 

The most natural approach in adapting finite-volume schemes to the MHD 
equations is to define the magnetic field component at the center of each cell, 
where the traditional hydrodynamical variables are also defined. One then 
takes advantage of decades of experience in the development of stable and 
accurate shock-capturing schemes. In this case, the solenoidality constraint 
V ■ B = has to be enforced using either a "divergence cleaning" step (see 
for example Brackbill and Barnes, 1980 and Ryu et al., 1998), or various 
reformulations of the MHD equations including additional divergence-waves 
(Powell et al., 1999) or divergence-damping terms (Dedner et al., 2002). A 
novel cell-centered MHD scheme has been recently developed by Crockett 
et al. (2005) that combines most of these ideas into one single algorithm. 

An alternative approach is to use the Constrained Transport (CT) algorithm 
for the induction equation, as suggested in the late 60's by Yee (1966), and later 
revisited by Evans and Hawley (1988). In this description, the magnetic field 
is defined at the cell faces, while other hydrodynamical variables are defined 
at the cell center. This is often called a "staggered mesh" discretization. As 
we will see in this paper, CT provides a natural expression of the induction 
equation in conservative form. Combining CT with the Godunov framework to 
design high-order, stable schemes is therefore a very attractive solution. This 
combined approach was first explored in the context of the MHD equations by 
Balsara and Spicer (1999). This method directly uses face-centered Godunov 
fluxes and averages these on the cell edges to estimate the Electro-Motive 
Force (EMF). Toth (2000) proposed an interesting cell-centered alternative 
to this scheme. More recently, Londrillo and Del Zanna (2000, 2004) have 
revisited the problem and shown that the proper way of defining the edge- 
centered EMF is to solve a 2D Riemann problem at the cell edges. They have 
applied this idea to design high-order, Runge-Kutta, ENO schemes. Finally, 
Gardiner and Stone (2005) have extended Balsara and Spicer scheme to design 
a more stable and more robust way of computing the EMF. 

The implementation of these various schemes within the Adaptive Mesh Re- 
finement framework is another challenging issue. It introduces two main new 

^ E-mail addresses: Romain.Teyssier@cea.fr (R.Teyssier), s.fromang@qmul.ac.uk 
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technical difficulties: ffist, proper fluxes and EMF corrections between dif- 
ferent levels of refinement must be accounted for. Second, when refining or 
de-refining cells, divergence-free preserving interpolation and prolongation op- 
erators must be designed. Both of these issues have recently been discussed in 
the framework of the CT algorithm by several authors (Balsara, 2001; Toth 
and Roe, 2002; Li and Li, 2004). 

The purpose of this article is to present a novel algorithm based on a high- 
order Godunov implementation of the CT algorithm within a tree-based Adap- 
tive Mesh Refinement (AMR) code called RAMSES (Teyssier, 2002). As op- 
posed to the grid-based (or patch-based) original AMR designed introduced 
by Berger and Oliger (1984) and Berger and Colella (1989), tree-based AMR 
trigger local grid refinements on a cell by cell basis. In this way, the grid fol- 
lows more closely the geometrical features of the computed flow, at the cost of 
a greater algorithm's complexity. Nevertheless, such tree-based AMR schemes 
have been implemented with success by various authors in the framework 
of astrophysics and fluid dynamics (Kravtsov et al., 1997; Khokhlov, 1998; 
Teyssier, 2002; Popinet, 2003) but not yet in the MHD context. On the other 
hand, patch-based AMR algorithms have been developed by several authors 
in recent years (Balsara, 2001; Kleimann et al., 2004; Powell et al., 1999; Sam- 
taney et al., 2004; Ziegler, 1999) and used for MHD applications. The main 
requirement that tree-based AMR usually place on the underlying solver is 
the compactness of the computational stencil: any high order scheme with a 
stencil extending to two points, or less, in each direction can easily be coupled 
to an "octree" data structure (Khokhlov, 1998). 

In this paper, our goal is to solve the induction equation using the MUSCL 
scheme, originally presented by van Leer (1977), and widely used in the litera- 
ture for the Euler equations. This very simple method is second order accurate 
in time and space and has a compact stencil: only 2 neighboring cells in each 
direction (and for each dimension) are necessary to update the central cell 
solution to the next time step. This compactness property is of particular im- 
portance for our tree based AMR approach. It is also useful for an efficient 
parallelization relying on domain decomposition. To our knowledge, this is the 
first implementation of the MUSCL scheme combined with the Constrained 
Transport algorithm that solves the induction equation. The key ingredient 
that ensures second order accuracy is the so-called "predictor step", in which 
the solution is first advanced by half a time step. We will consider a few differ- 
ent computational strategies for this predictor step and discuss their respective 
merits. Finally, we will present our overall tree-based AMR scheme. 

This paper is limited to the induction equation. We intend to apply the same 
approach to the full MHD equations in a future paper. Nevertheless, it is 
interesting to determine if such a numerical approach can be applied to kine- 
matic dynamo problems, for which the induction equation alone applies. The 
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induction equation is linear, but it can yield remarkably rich magnetic insta- 
bilities corresponding to exponential field growth and referred to as "dynamo 
instabilities" . The description of these instabilities, and the conditions under 
which they occur, constitute an active field of research, with important conse- 
quences in astrophysics and in geophysics, since they account for the origin of 
magnetic fields in the Earth, planets, stars and even galaxies. We will restrict 
our attention here to well known dynamo flows and use them to investigate 
the numerical properties of our scheme. 



An important problem in dynamo theory is related to a subclass of dynamo 
flows, known as "fast dynamos" which yield exponential field growth with 
finite growth rates in the limit of vanishing resistivity. This is of particular 
importance for astrophysical applications. Fast dynamos, when investigated 
with small, but finite, resistivity yield eigenmodes that are very localized in 
space, and are therefore ideal candidates for an investigation using the AMR 
scheme. 



Dynamo problems have traditionally been studied using spectral methods 
(Galloway and Frisch, 1986; Christensen et al., 2001). Some recent models 
have been produced using finite differences (Archontis et al., 2003), finite vol- 
umes (Harder and Hansen, 2005) or finite elements (Matsui and Okuda, 2005). 
However, all of these methods rely on explicit physical diffusion to ensure nu- 
merical stability. The interest of using CT within the Godunov framework 
together with an AMR approach is twofold. First, fast dynamo modes have a 
very localized spatial structure (scaling as Rm~^^'^ where Rm is the magnetic 
Reynolds number). Adapting the computational grid to the typical geome- 
try of these modes therefore appears as a very natural strategy to minimize 
computational cost. Second, the Godunov methodology, using the CT scheme, 
introduces the minimal amount of numerical dissipation needed to ensure sta- 
bility. This is an important property when using an AMR approach, for which 
cells of very different sizes coexist. This last property of the scheme is then 
mandatory to allow the use of a coarse grid in regions barely affected by the 
physical diffusion. 



We will present several tests that demonstrate the efficiency of our tree-based 
AMR Godunov CT scheme for solving complex dynamo problems: we will first 
reproduce a simple advection problem of a magnetic loop and then validate 
the approach on two well known dynamo fiows: the Ponomarenko dynamo and 
a fast ABC dynamo. 
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2 Constrained Transport in Two Space Dimensions 



In this section, we briefly review the design of stable numerical schemes for 
hyperbolic systems of conservation laws in two space dimensions using the 
Godunov approach. Following Londrillo and Del Zanna (2000), such systems 
are called here "Euler systems" , as opposed to the "induction system" we will 
consider later. 



2.1 First Order Godunov Scheme for Euler systems 



We flrst examine the problem in one space dimension. The following Euler 
system, 

— + V-F(U)=0, (1) 



can be written in integral form by defining finite control volume elements in 
space and time, where we define a cell by Vi = [x , i , a^.^i] and a time interval 

by At = t^^^ — t^. The conservative system writes for each cell Vi 

At ' 1 



Note that this integral form is exact for the corresponding Euler system. The 
averaged, cell-centered state is defined by 

X 1 

J U{xX)dx, (3) 



while the averaged, time-centered intercell flux is deflned by 



The Godunov method states that the intercell fiux is computed using the 
solution of a Riemann problem with left and right states given by the left and 
right averaged states 



u:,_{x/t)^Rp \{u):,{u)i 



i+ 



i+1 



(5) 
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This approach, called "first order Godunov scheme" , assumes that the solution 
inside cell Vi is piecewise constant. Taking advantage of the self-similarity of 
the Riemann solution for initially piecewise constant states, one can simplify 
further the time-average of the fiux and obtain 

f"^,^ = f{u;ao)) . (6) 



Note that again the time evolution of the average state over one time step 
is exact. Numerical approximations arise when one assumes at the next time 
step that the new solution inside cell Vi is also piecewise constant and equal 
to the new averaged state. 

We now extend the previous method to Euler systems in 2 space dimensions. 
The conservative system can also be written in the following unsplit formula- 
tion 

- m + ^ i F^'i - F^'i ] + ^ (g-^^ - g^'^K] = 0,(7) 



where the average state is now defined over a 2 dimensional cell Vij, and 
intercell fluxes are now time averaged fiuxes integrated over the line separating 
neighboring cells 



y 1 



!/ 1 

'-2 



2 t" X 1 

'"2 



At this point, the integral form is still exact. The generalization of the ID 
Godunov scheme to multidimensional problems now relies on solving two di- 
mensional Riemann problems at each corner, defined by four initially piecewise 
constant states 



U* 



{x/t,y/t) = RP f([/)^ , {U)l,_^ {u)r^, {U)l 



'i+lj+lj 



(10) 



The fundamental difference with the ID case is that we now need to average 
the complete solutions of 2 adjacent Riemann solutions over the entire trans- 
verse fine segment, where fiuxes are defined. These space- averaged fiuxes are 
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not functions of a unique self-similar variable anymore, but depend explicitly 
on time. Building such a numerical scheme is barely possible for simple scalar 
linear advection problem and far too complex to implement for non-linear 
systems. 

The traditional approach is to approximate the true solution using a predictor- 
corrector scheme. This is also the key ingredient of any high-order scheme, 
where the self-similarity of the Riemann problem breaks down, even in one 
space dimension, due to the underlying piecewise linear or parabolic repre- 
sentation of the data. The idea is to compute a predicted state at time level 
use this intermediate state as an input state for the two final ID 
Riemann solvers. 

We list here 3 classical methods to implement this predictor step 

• Godunov method: no predictor step is performed. This greatly simphfies 
the method, which now relies on one Riemann solver in each direction. The 
prize to pay is a somewhat restrictive Courant stability condition: {u/Ax-\- 
v/Ay)At < 1, where u and v are the maximum wave speed in each direction. 

• Runge-Kutta method: the predictor step is performed using the 2D Go- 
dunov method with half the time step. The resulting intermediate states 
are then used to compute the fluxes for the final conservative update. The 
Courant condition is the same as for the Godunov method, but one has to 
perform 2 Riemann solvers per cell in each direction (4 in total). 

• Corner Transport Upwind method: predicted states for a given Rie- 
mann problem are computed with a ID update in the transverse direction 
only, for the time interval At/2. This scheme was first proposed by Colella 
(1990). It allows up to a factor of two larger time steps than the two previ- 
ous schemes, since the Courant condition is now max{u/ Ax,v/ Ay)At < 1, 
but 2 Riemann solvers per cell in each direction (4 in total) are still needed. 

All three methods are directionally unsplit, first order approximations (in 
space) of the underlying Euler system. 



2.2 First Order Godunov Scheme for the Induction Equation 



The magnetic field evolution in the MHD approximation is governed by the 
induction equation which neglects free charge density and displacement cur- 
rents. It is written in conservative form as 

-^ = VxE + 77AB, (11) 
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where the EMF E is given by 

E = vxB, (12) 



and 7] is the magnetic diffusivity. The magnetic field also satisfies the diver- 
gence free constraint 

V-B = 0. (13) 



It is usually more convenient to consider (11) in non-dimensional form by 
introducing a typical lengthscale C and a typical timescale T — C/U where U 
is some norm of the velocity (usually based on the maximal value over space 
and time). The resulting non-dimensional equation is 

= V X (v X B) + Rm-^ AB, (14) 



where Rm = {UC)/rj while t = t/T and v = ^^/U are respectively the non- 
dimensionnal time and velocities and the spatial derivative are taken with 
respect to normalized distances. 

The EMF E is here the analog of the flux function for Eulcr systems. Wc now 
restrict our attention to 2D dimensional flows'^ , for which only one component 
of the EMF, say E^, is sufficient. 

Following the Godunov approach, we write the 2D induction equation in in- 
tegral form over a finite control volume in space and time. For the com- 
ponent of the magnetic field, we define a finite surface element 5'i+i/2j = 
[yj_i/2,%+i/2] at position Xi+1/2 

=(5.);i +^f(EX;i^^i- 1) . (15) 



For the By component, we define a finite surface element S'jj+1/2 = [a;j_i/2, 2:^+1/2] 
at position yi+i/2- The induction equation in integral form has a similar rep- 
resentation 

{By}:^li = {By):^^i - ^ ({E^c'i .1 - {Es:i ^ . (16) 

hJ+2 *'.J + 2 l^X \ 1+2^3 + 2 *~2'-' + 2/ 

^ The one dimensional induction equation, with =constant, is equivalent to a 
Euler system, for which the standard methodology applies without modification. 
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Note that this integral form in space and time is exact. The average, surface 
centered, magnetic states are defined as the average magnetic field components 
on their corresponding control surfaces 



y 1 



<-^'>Hi> = ij / BJ^.i'y-ndy, (17) 



Ay 

1 

'-2 



^^^.4 = 2^ I By(-^y.4^ndx. (18) 



X 1 

'2 



2.2.1 2D Riemann Problem 



The time centered EMF results from a time average at the corner points 

{E,)''y f E,{x i,y i,t)dt. (19) 



Let us now apply the Godunov method to the 2D induction equation. Upon 
noticing that our initial conditions are given by four piecewise constant states 
around each corner points, we can use the self-similar solution of the 2D Rie- 
mann problem at the corner point. 



U* 



{x/t,y/t) = RP [([/)^ , {U)l,^^ {U)r^, {U)l 



(20) 



and time integration vanishes in equation (19) 

»+2'J + 2 *+2'-^'^2 



The Godunov method, applied to the induction equation in 2D, shares this 
interesting property with the Godunov method applied to ID Euler system. 
The self-similarity of the flux function was lost for 2D Euler systems. The 
self-similarity of the EMF function is still valid for the 2D induction equation, 
provided our initial conditions are described by piecewise constant states. We 
will see in the next section, that this is unfortunately not true in the general 
case, even at lowest order. 

As noticed by Londrillo and Del Zanna (2000), the 2D Riemann problem is 
the key ingredient for solving the induction equation with a stable (upwind) 
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Fig. 1. The 2D Riemann problem in the x-y plane to compute the EMF in the z 
direction at edge {i + ^, j + The face-centered magnetic fields are shown as 
vertical and horizontal arrows. The velocity field is shown as the dashed arrow. 

scheme. The 4 initial states (with 2 magnetic field components per state) need 
to satisfy the V • B = property. should therefore be the same for the 
two top states, and for the two bottom states, while By should be the same 
for the two left states, and for the two right states (see Fig. 1). This condition 
is naturally satisfied as long as magnetic field is defined as a surface-average, 
see (17) and (18). 

In the general MHD case, designing 2D Riemann solvers (even approximate 
ones) is a very ambitious task. For the kinematic induction case, the solution is 

however remarkably simple, since the solution is nothing else but the upwind 
state. The edge-centered EMF can therefore be written in the following closed 
form 



\^z) 1 1 — U V 

-H -1^ ^ + ^ (22) 

where u and v are respectively the x and y components of the flow velocity 
V = (m, V, w) computed at the center of the edge {i + |, J + |). This last equa- 
tion is familiar in the framework of upwind finite-volume schemes. It can be 
decomposed into two contributions. The first line is the EMF computed using 
the average magnetic fields at the cell corners: this EMF is a second-order 
in space. The resulting scheme (retaining this term only) would have been 
unconditionally unstable, if it was not for the second term, the contribution 
of the upwinding. It is equivalent to a 2D numerical diffusivity, with direc- 
tional diffusivity coefficients given by r]^ = \u\ Aa;/2 and rjy = \v\ Ay/2. This 
(relatively large) resistivity introduces the minimal but necessary amount of 
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numerical diffusion for the scheme to remain stable. 



2.2.2 Constrained Transport as a Finite Surface Approximation 

This straightforward extension of the Godunov methodology has lead us to 
the well known "Constrained Transport" (CT) scheme, that was designed a 
long time ago for the MHD equations by Yee (1966). The key property of the 
CT scheme is that one can also write the V • B = constraint in integral form 

as 



This integral form is exact. Moreover, if it is satisfied by our initial data, the 
integral forms in (15) and (16) ensure that it will be satisfied at all iterations 
during the numerical integration. Using equation (23), and assuming that 
formally Ax ^ 0, we show that the following property holds: 

Remcirk 1 {B^)^. (x) is a continuous function of coordinate x, 

and, symmetrically, assuming that formally Ay — > 0, we have: 

Remark 2 {By)^ (y) is a continuous function of coordinate y, 

This means that (Sa;)"^^^2j can be considered as piecewise constant in the y 
direction, but has to be considered as piecewise linear in the x direction. This 
constitutes our lowest order approximation of the magnetic field. Symmetri- 
cally, to lowest order, (-Bj,)"j_|_i/2 ^^"^ considered as piecewise constant in the 
X direction, but has to he considered as piecewise linear in the y direction. ^ 

This last property provides a fundamental difference between the induction 
equations and Euler systems. It is due to the divergence free constraint, ex- 
pressed in integral form on a staggered magnetic field representation. One 
consequence of this property is that our initial state for the 2D Riemann 
problem cannot be piecewise constant anymore, but instead piecewise linear. 
We therefore loose the property of self-similarity for the Riemann solution at 
corner points, and cannot perform an exact time integration to compute the 
time average EMF. We now have to rely on approximations. Following the 
strategies developed in section 2.1, we approximate the time averaged EMF 
using various predictor-corrector schemes. 



^ Let us stress that for ideal MHD, a jump perpendicular to the fieldline is allowed. 




Ax 



Ay 



0. 



(23) 
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2.3 The Predictor step 
2.3.1 Godunov Scheme 

The first possibility is to drop the predictor step and solve the Riemann prob- 
lem defined at time f^. Using (15) and (16), together with the EMF computed 
from (22), we obtain the Godunov scheme for the induction equation. In the 
simple case of a constant velocity field with u > Q and v > Q (the pure 
advection case), we can write the overall scheme as 



Using the V • B = constraint at time in integral form (23), we further 
simplify the scheme to obtain 



We can therefore conclude: 

Proposition 1 For the advection case, if the initial data satisfy the integral 
form of the solenoidality constraint, the Godunov method for the induction 
equation is identical to the Godunov method for the advection equation on the 
staggered grid. 

This rather simple point is actually quite important, since it proves that CT 
has advection properties quite similar (in this case identical) to traditional 
finite-volume methods. The Godunov scheme for the induction equation has 
a compact stencil. It is however of mere theoretical interest, since, as we will 
see in the next section, it is not the first order limit of higher order Godunov 
implementations of the induction equation. 

2.3.2 Runge-Kutta Scheme 

As discussed above, the V ■ B = constraint, and the loss of self-similarity in 
the Riemann solution, pushes towards using a predictor step in designing our 
first order scheme. The most natural approach is the Runge-Kutta scheme, 
for which the solution is advanced first to the intermediate time coordinate 




(24) 




(25) 
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Fig. 2. Stencils of our various schemes for the induction equation: Runge-Kutta 
scheme (left plot), U-MUSCL scheme (middle plot) and C-MUSCL scheme (right 
plot). The flux being computed is indicated by a bold face and arrow. For the pur- 
pose of this example, the velocity field is pointing in the upper right direction (u > 
and f > 0). The first order stencil in space (second order in time) is represented 
with black arrows. Additional components required for the second order stencils in 
time and space are shown with white arrows. The shaded region indicates cells that 
are available in a tree-based AMR implementation. Only the two right schemes have 
stencils compact enough for such an implementation. 

using the (previously described) Godunov scheme with time step At/ 2. 
These predicted states are then used to define the 4 initial states for the 2D 
Riemann problem. The resulting EMF is used to advance the solution from 
time to the next time coordinate t""*"^ with time step At. A similar, 2 
step, Runge-Kutta method for the induction equation is used for example in 
Londrillo and Del Zanna (2000) and Londrillo and Del Zanna (2004) to solve 
the MHD equations. 



Using similar arguments as in the previous section, it is easy to show that, for 
a uniform velocity field, since the predicted magnetic field satisfies the inte- 
gral form of the solenoidality constraint, the corrector step for the induction 
equation is identical to the predictor step for the advection equation. As we 
have shown in the last section, this property also holds for the predictor step, 
we therefore obtain a second important result: 

Proposition 2 For a uniform velocity field, if the initial data satisfy the in- 
tegral form of the solenoidality constraint, the Runge-Kutta method for the 
induction equation is identical to the Runge-Kutta method for the advection 
equation on the staggered grid. 

We will show later that it is also possible to design higher order schemes 
for this algorithm. This scheme has two nice properties: it is second order in 
time (while still first order in space), and the predicted magnetic field satisfies 
exactly V • B"+^/^ = 0. There are also issues associated with it, especially in 
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the AMR framework. It can be easily shown (see Fig. 2) that the stencil is 
not compact enough for a tree-based AMR: 3 ghost cells are needed in each 
direction (resp. 2) for the second order (resp. first order) scheme. We will see 
in the test section that it is also shghtly more diffusive than the other schemes 
we will describe in the following sections. The Courant stability condition is 
also rather restrictive 

2.3.3 Upwmd-MUSCL Scheme 

When deriving the MUSCL scheme for Euler systems, van Leer (1977) noticed 
that it was not necessary for the predictor step to be strictly conservative. A 
conservative update was however mandatory for the corrector step. Similarly, 
for the induction equation, it is a priori not necessary for the predictor step 
to satisfy the solenoidality constraint. It is however mandatory for the initial 
and final data. Instead of computing one EMF at each cell corner, using a 
2D Riemann solver, we now propose to compute for the predictor step only 4 
EMFs at each cell corner, corresponding to each input magnetic field. 

These EMFs are defined as (^,)J;i/2,j+i/2' (^^) Ji/2j+i/2' {Ez)f+i/2,j+i/2 and 
j+i/2j+i/2' where each upper index corresponds to the "left", "right", 
"bottom" and "top" face, respectively. Each EMF is speciahzed to its cor- 
responding face-centered magnetic field component. One EMF per face is al- 
lowed, in order to satisfy the continuity constraint: we need to solve a ID 
Riemann problem in the perpendicular direction. The Riemann solution is 
here the "upwind" state. The "bottom" and "top" EMF for the predictor step 
are therefore 




Similarly the "left" and "right" EMF are 
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+ l''l((B.^,„-(B.).,l,)/2- (28) 
The predictor step for the x component of the magnetic field becomes 

{B^ri/' = (5.);, + ((£;.)f i 1 - 1 i ) , (29) 

and for the y component we have 

(5.)"+^' = (S,)" i-4^f(E,)^i iV (30) 

\ y^^j+l \ 2AX V ''+2'^+2 ^ ^'-2'J + 2/ ^ ^ 

To complete this scheme, the corrector step is performed using a final 2D Rie- 
mann solver to compute the time-centered EMF (22) and a final conservative 
update of each magnetic field component (15) and (16). 

Let us now examine the property of the Upwind-MUSCL scheme in the case 
of a uniform velocity field. We can assume, without loss of generality, that 
u > and v > 0. In this case, the predicted state can be written in a more 
compact form 

which is equivalent, using (23), to 

= (S,)" 1 - ( (S,)" 1 - (5,)" 1 ) . (32) 

Similar expressions can be derived for {By)'^^^y2- Inserting these predicted 
values into (22) and (15), we get, after some tedious manipulations, the final 
updated solution 

{B^r^l . = {B^l^ (1 - a) (1 - Cy) + {B^y; 1 a. (1 - Cy) 

+ (1 - + ^^-)r-i j-1 ^^^y ' (33) 

where the following definitions have been used = uAt/Ax and Cy — 
vAt/Ay. One can recognize here the Corner Transport Upwind (CTU) ad- 
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vection scheme presented in Colella (1990), for which the Courant stabihty 
condition is 



u 



max 



Ax' Ay 



At<l. 



(34) 



We therefore conclude: 

Proposition 3 For a uniform velocity field, if the initial data satisfy the in- 
tegral form of the solenoidality constraint, the Upwind-MUSCL Scheme for 
the induction equation is identical to Colella 's first order CTU scheme for the 
advection equation on the staggered grid. 

It is apparent in (33) that the stencil of this MUSCL scheme is more compact 
that it is for the Runge-Kutta scheme (see also Fig. 2). Since our goal is here 
to develop an AMR code for the induction equation, this is a very attractive 
solution. The predictor step is performed using upwinding in the normal di- 
rection. As for Colella's CTU scheme, the Courant stabihty condition is very 
efficient. We now explore one last possibility for our MUSCL predictor step. 



2.3.4 Conservative-MUSCL Scheme 

The last scheme was designed in dropping the solenoidality constraint for the 
predictor step. We propose in this section to drop the upwinding in the EMF 
computation for the predictor step, which now becomes 

{E.)ll^^^l = u '^^^ ^ - . '^^^ ^ . (35) 



Since wc now have a single EMF per cell corner, the predicted magnetic field 
satisfies by construction V ■ B"^^/^ = 0. The corrector step is the same as 
for all 3 methods. Here again, we would like to examine the property of the 
scheme for the case of uniform advection. Because in this case V • B"'^^''^ = 0, 
the corrector step is identical to the corrector step for the Godunov advection 
scheme on the staggered grid. The predictor step, on the other hand, can 
be written as the Forward Euler scheme for the advection equation on the 
staggered grid. When combined together, we obtain a new first order advection 
scheme for which the Courant stability condition is the same as for the Runge- 
Kutta scheme. For this new scheme to be monotone, however, the time step 
has to satisfy the following more restrictive condition 

\Ax Ay -^2 + 1 ^ ' 
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Proposition 4 For a uniform velocity field, if the initial data satisfy the inte- 
gral form of the solenoidality constraint, the Conservative-MUSCL Scheme for 
the induction equation is identical to a new, consistent and stable first order 
scheme for the advection equation on the staggered grid. 

At the expense of a more restrictive constraint on tlie time step, we have 
obtain a new scheme which is conservative for the predicted step, in the sense 
that the predicted magnetic field satisfies the solenoidality constraint. 



2.4 High Order Schemes 



Extensions of the three above schemes (Runge-Kutta, U-MUSCL and C- 
MUSCL) to second order are based on a piecewise linear reconstruction of 
each magnetic field component, using "magnetic flux conserving" interpola- 
tion at each cell interface. Following the MUSCL approach, one can compute 
corner (or edge) centered interpolated quantities, using a Taylor expansion 
both in time and space as follows, for 



\n+l/2,B 
\-^x/ 11 

/p \n+l/2,T 
\^x/ 1 .1 
2+ 2 iJ 2 



^+2.J 



^dB, 
K dt 



\ — + 



At 



'ae^r Ay 



dy 



^. (37) 



and for B,, 



I rj \ n+l/2,L 

\-^y/ _. , 1 . , 1 




-IT ■ (38) 



In this way, second-order, edge-centered components of the magnetic field 
can be used in the 2D Riemann solver to compute the EMF and update 
the solution to time Our three different schemes differ in the way they 
implement the terms dB^/dt and dBy/dt. 

Let us stress that to recover second order accuracy in space, one needs to 
perform a predictor step which is also second order accurate in space. For the 
C-MUSCL scheme, this is already the case if one uses exactly the predictor step 
presented in the last section. For both the Runge-Kutta and the U-MUSCL 
schemes, however, one needs to use a linear reconstruction of each magnetic 
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field component and compute the EMF for the predictor step. This is done 
using the following equations 



These edge-centered components are then used to compute the EMF, using 
(22) for the Runge-Kutta method, or (27) and (28) for the U-MUSCL scheme. 
As usually done in higher order finite volume schemes, spatial derivatives are 
approximated using slope limiters, in order to obtain positivity preserving, non 
oscillatory solutions. For that purpose we use a standard slope limiter (used in 
many fluid dynamics codes), the Monotonized Central Limiter, which is given 
by 

Far from discontinuities, this slope reduces to Fromm's finite difference ap- 
proximation of the spatial derivative. In this case, one can show that, for 
a uniform velocity field, all 3 schemes are again strictly equivalent to their 
second order parent scheme for the advection equation on the staggered grid. 

In non smooth parts of the flow, however, this is no longer true. Slope lim- 
iting destroys the strict equivalence between the induction schemes and their 
advection counterparts. One must also be aware that traditional slope lim- 
iters, such as the one we use here, are designed for the advection equation 
in finite-volume schemes. The monotonicity of the solution for the induction 
equation is therefore not guaranteed. Deriving slope limiters for the induction 
equation is beyond the scope of this paper. We have to rely on the numerical 
tests performed in the test section to assess the non oscillatory properties of 
our schemes. 

It is also apparent in (41) that for both Runge-Kutta and U-MUSCL schemes, 
the computational stencil increases by one cell in each direction, compared 
to the first order scheme (see Fig. 2). The second order U-MUSCL and the 



18 



C-MUSCL schemes are therefore both compact enough for our AMR imple- 
mentation, while the second order Runge-Kutta scheme is not. 

2.5 Conclusion 

We have derived in this section three numerical schemes for the solution of the 
induction equation using the CT algorithm in two-dimensions. All of them are 
second order in space and time. We have called these schemes Runge-Kutta, 
U-MUSCL and C-MUSCL. Only the last two have compact computational 
stencils, which makes them suitable for our tree-based AMR implementation. 
More interestingly, we have proven that, in case of a uniform velocity field, the 
U-MUSCL scheme is strictly identical to Colella's Corner Transport Upwind 
scheme for the advection equation on the staggered grid. For the C-MUSCL 
scheme, we have shown that it is strictly identical to another wcll-behaved ad- 
vection scheme, with however a more restrictive stability condition on the time 
step. This shows that CT, when properly derived within Godunov's frame- 
work, has advection properties similar to traditional finite-volume schemes. 



3 A Constrained Transport AMR Scheme in three Dimensions 

In this section, we describe our MUSCL-type schemes for the induction equa- 
tion in three space dimensions. It is mostly a straightforward generalization 
of the previous 2D schemes, we will however repeat each step of the algorithm 
in order to summarize our method, and introduce the discussion of the AMR 
implementation. 

3.1 Definitions 

Let us generalize the schemes discussed in 2D in section 2 to 3D problems. 
The three magnetic field components are discretized on a staggered grid using 
a finite-surface representation 

= / / 5.(^m/2,y,.,ndz/d., (42) 

2/i-l/2 Zi-1/2 

^i+1/2 ^1+1/2 

= / / By{x,y,^,/2,zX)dxdz, (43) 
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^""^K.^^^^^ I I B.ix,y,z,,y,,nd.dz. (44) 

These three conservative variables satisfy the divergence-free constraint in 
integral form 

5.^ Conservative update 

The magnetic field components are updated from time to time using 
the induction equation in integral form, which becomes (for B^) 

see (15) for comparison. 

Similar expressions can be derived for By and B^. Here, E'^;, Ey and E'^ are 
time-averaged EMFs defined at each cell edges. 

3.3 2D Riemann Solver 

Each of these EMFs components arc obtained as the solution of a 2D Riemann 
problem, defined by 4 initial states surrounding the corresponding edge. The 
upwind solution of this 2D Riemann problem for is given by 

y «,:;+2'«^+2 »j+2''^+2/ 
20 



(45) 



Where the magnetic field components, labeled n+1/2, R; n+1/2, L; n+1/2, T 
and n+1/2, B are the time-centered predicted states interpolated at cell edges. 
Similar expressions for Ey and can be deduced by permutations. 



3.4 Predictor Step 



The predicted states of the magnetic field are obtained through a Taylor ex- 
pansion in time and space. For B^, this translates into 



, ,„+i/2_B /ID \n ( dB^y^ At ( dB^y^ Ay 



2' 

n 



, n+1/2, _B _ / n \n . ( dBx\^ At . / dB^\^ Az 



Similar expressions can be written for By and B^ . The spatial derivatives are 
computed in each direction using the slope limiter function (41). Our three 
schemes differ only in the way the time derivative is estimated in the above 
expansion. 



3.4-1 Runge-Kutta Scheme 

The Runge-Kutta predictor step is equivalent to the corrector step, except 
for the time derivative in (48). We use spatial derivatives to define edge- 
centered magnetic field components and the 2D Riemann solver to define the 
edge-centered EMF components. This unique EMF vector, defined at time t", 
is finally used in the conservative formula (46) to obtain a finite difference 
approximation of the time derivative in (48). For a uniform velocity field, 
the first order scheme is again identical to the Runge-Kutta scheme for the 
advection equation on the staggered grid. For the second order scheme, this 
is only true in smooth regions of the solution. 
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3.4-2 U-MUSCL Scheme 



For the U-MUSCL scheme, the EMF used to compute the predicted states is 
not uniquely defined at each edge anymore, so that the predicted magnetic 
field does not satisfy the divergence-free constraint. In fact, we compute at 
each cell edge 4 EMF components, specialized to each face-centered magnetic 
field component. By solving a ID Riemann problem at each faces, we perform 
a proper upwinding in the normal direction. The input states of these ID 
Riemann problem are reconstructed magnetic field components at cell edges 
using slope limiters. Note that for a uniform velocity field, this first order 
scheme is not equivalent anymore to the CTU scheme in 3D. 

3.4.3 C-MUSCL Scheme 

Like the Runge-Kutta method, the C-MUSCL scheme involves one single EMF 
vector to compute the time-derivative in the Taylor expansion, therefore pre- 
serving the solenoidal property on the predicted step. This EMF is computed 
using the average of the face-centered magnetic field components, as in (35). 
It does not involve any limited slope computations, but still retains second 
order accuracy in space. As explained in the previous section, the cost is a 
more restrictive time-step stability condition. For a uniform velocity field the 
scheme is identical to the new advection scheme on the 3D staggered grid 
discussed in section 2.3.4. 

3.4.4 Merits of the Various Schemes 

We compare, in this section, the different advantages and drawbacks of each 
of the above described methods. The corrector step is the same for each cases. 

The Runge-Kutta scheme is the most natural scheme to write. However, it will 
prove to be very expensive for MHD, since it requires a 2D Riemann solver 
in the predictor step. Moreover, it has a restrictive Courant condition and its 
stencil is too large to be implemented in the AMR implementation, which is 
not the case of the two other schemes. 

The U-MUSCL scheme has better stability properties, the time step is less 
restrictive. It is also expected to be more efiicient in MHD applications, since 
one ID Riemann problem only is required in the predictor step. Note however 
that its rigorous 3D extension is problematic and requires further investigation. 

Unlike the U-MUSCL scheme, for which the non-conservation of the solenoidal- 
ity condition in the predictor step may cause problems in some cases, the C- 
MUSCL scheme is conservative. No Riemann solver is needed in the predictor 
step, which should make it very efficient for MHD (Fromang et al. (2006)). 
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But these advantages are obtained at the cost of a smaller timestep than the 
U-MUSCL scheme. 



3.5 AMR Implementation 

We have included both of the compact schemes (U-MUSCL and C-MUSCL) 
in the RAMSES code. It is a tree-based AMR code originally designed for 
astrophysical fluid dynamics (Teyssier, 2002). The data structure is a "Fully 
Threaded Tree" (Khokhlov, 1998). The grid is divided into groups of 8 cells, 
called "octs" , that share the same parent cell. Each oct has access to its parent 
cell address in memory, but also to neighboring parent cells. When a cell is 
refined, it is called a "split" cell, while in the opposite case, it is called a 
"leaf" cell. The computational domain is always defined as the unit cube, 
which corresponds in our terminology to the first level of refinement in the 
hierarchy £ = 1. The grid is then recursively refined up to the minimum level 
of refinement £mmi in order to build the coarse grid. This coarse grid is the 
base Cartesian grid, covering the whole computational domain, from which 
adaptive refinement can proceed. This base grid is eventually refined further 
up to some maximum level of refinement Imaxi according to some user defined 
refinement criterion. 

When Irnax = ^min, the Computational grid is a traditional Cartesian grid, 
for which the previous induction schemes apply without any modification. 
When refined cells are created, however, some issues specific to AMR must be 
addressed. 



3.5.1 Divergence-free Prolongation Operator 

When a cell is refined, eight new cells (i.e. a new "oct") are created for which 
new magnetic field components are needed. More precisely, each of the six faces 
of the parent cell are split into 4 new fine faces. Three new faces, at the center 
of the parent cell, are also spht into four new children faces. The resulting 
magnetic field components, fine or coarse, need to satisfy the divergence-free 
constraint in integral form. 

This critical step, usually called in the multigrid terminology the Prolongation 
Operator, has been solved by Balsara (2001) and Toth and Roe (2002) in the 
CT framework. We recommend both of these articles for a detailed description 
of the method. The idea is to used slope limiters to interpolate the magnetic 
field component inside each parent face, in a flux-conserving way, and then to 
use a 3D reconstruction, which is divergence-free in a local sense inside the 
whole cell volume, in order to compute the new magnetic fleld components 
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for each central children faces. In our case, the same slope limiter as in the 
Godunov scheme (41) has been used. 

This prolongation operator is used to estimate the magnetic field in newly 
refined cells, but also to define a temporary "buffer zone" , two "ghost cells" 
wide, that set the proper boundary for fine cells at a coarse- fine level boundary. 
This is the main reason why a compact stencil is needed for the underlying 
Godunov scheme. 

3.5.2 Magnetic Flux Corrections 

The other important step is to define the reverse operation, when a split 
cell is de-refined, and becomes a leaf cell again. This operation is usually 
called the Restriction Operator in the multigrid terminology. The solenoidal- 
ity constraint needs again to be satisfied, which translates into conserving 
the magnetic fiux. The magnetic field component in the coarse face is just the 
arithmetic average of the 4 fine face values. This is reminiscent of the "flux cor- 
rection step" of AMR implementations for Euler systems (Berger and Ohger, 
1984; Berger and ColeUa, 1989; Teyssier, 2002). 

3.5.3 EMF Corrections 

The "EMF correction step" is more specific to the induction equation. For a 
coarse face which is adjacent, in any direction, to a refined face, the coarse 
EMF in the conservative update of the solution needs to be replaced by the 
arithmetic average of the two fine EMF vectors. This guarantees that the 
magnetic field remains divergence-free, even at coarse-fine boundaries. 

3.6 Physical resistivity 

We have now completely described our AMR implementation for the induction 
equation. It can be used as such, without explicitly including physical resis- 
tivity, to investigate fast-dynamo action associated with a given flow. The 
resulting integration is stable and produce an exponentially growing field very 
similar to what we expect in dynamo theory. However, resistivity (and thus 
reconnection) , which is necessary to identify a growing eigenmode, is solely 
due to the underlying numerical scheme. This numerical resistivity is usually 
non-uniform in time and space, anisotropic and non-linear. The mathematical 
properties of the resulting eigenmodes are unclear, and the results usually de- 
pend on the mesh resolution. Instead, we have chosen to explicitly introduce 
a physical resistivity in the induction equation, see (14), in order to allow a 
proper identification of the eigenmode. 
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The amplitude of the resistive term is here controlled by the inverse of the mag- 
netic Reynolds number Rm = UL/rj. We shall concentrate on large magnetic 
Reynolds numbers (i.e. the fast dynamo limit). It may, at first, seem strange 
to introduce this term when the Godunov approach has precisely been intro- 
duced to ensure numerical stability and reduced numerical diffusion. In fact, 
because of the very nature of the fast dynamo solution, the effect of physical 
resistivity will be limited to very localized regions. Its effect will therefore be 
limited to the very fine AMR cells and the stabilizing property of the Godunov 
approach will be essential for the coarser cells. 

Physical diffusivity is introduced in our scheme using the operator splitting 
technique. After the induction equation has been advanced to the next time 
coordinate with solution B*, we solve for the diffusive source term, using 
the following equation 

■pn+l TD* 

= ,yVxj"+^ where = V x B"+\ (49) 



where j is the current. It is defined at cell edges. For example, the finite 
difference approximation for {jy and jz are not shown) is written as 



-{B. 



{By 



-{By) 



X). .,1 1 



n+1 



At/ 



(50) 



Considering the current as the analog of the EMF, all the ingredients of the 
previous sections can be applied to design a conservative AMR implementa- 
tion to solve for the diffusion source term. We use for that purpose a fully 
implicit time discretization, in order for the time step to be limited only by 
the induction scheme Courant stability condition. The resulting linear sys- 
tem is solved iteratively using the Jacobi method. Note that in the problems 
we address in this paper, only a few iterations were necessary to reach 10~^ 
accuracy. 



4 Tests and Application to Kinematic Dynamos 



In this section, we test our various schemes using the advection of a magnetic 
field loop in 2D. We conclude that the three Godunov schemes we described 

for the induction equation have very good and similar performances. The U- 
MUSCL scheme seems to be slightly better than the other two. We also test the 
AMR implementation, showing that the results are almost indistinguishable 
from the reference Cartesian run. We will then use this code to compute the 
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Runge-Kutta 1st order Runge-Kutta 2nd order 




Fig. 3. Magnetic loop advection test for a Cartesian grid with Hx = 128 and Uy = 64: 
each panel shows a gray-scale image of the magnetic energy (B^ + By) at time t = 2. 

The scheme used to compute each image is provided in the title of each panel. Sec- 
ond-order schemes give very similar results, while the first order U-MUSCL scheme 
performs slightly better than the two other first order schemes. 

evolution of two well-studied dynamo flows: the Ponomarenko dynamo and 
the ABC flow. This will serve as a final integrated test of our scheme. 

4-1 Magnetic Loop Advection 

Let us first focus our attention on a simple test of pure advection which was 
recently proposed by Gardiner and Stone (2005) to investigate the advection 
properties of their CT scheme. It consists in the advection of a magnetic field 
loop with a uniform velocity field. It is of particular relevance in our case, since 
we are dealing with kinematic induction problems. The computational domain 
is defined by — 1 < x < 1 and —0.5 < y < 0.5. The boundary conditions are 
periodic. The flow velocity is set to = 2, v = 1 and w — 0. 

The initial magnetic field is such that = 0, while and By are defined 
using the z-component of the potential vector A (with B = V x A), as an 
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Fig. 4. Magnetic energy as a function of time for the field loop advection test. The 
upper solid line is the solution for perfect advection. The lower lines are for the first 
order schemes: Runge-Kutta (dotted line), C-MUSCL (dashed line) and U-MUSCL 
(solid line). Runge-Kutta and C-MUSCL results are indistinguishable in this case. 
The 3 intermediate lines correspond to second order schemes and use the same line 
convention. The dot-dashed lines is the AMR result obtained with U-MUSCL and 
using irnin = 3 and imax = 9- 

axisymmetric function of the form 
R — r for r < R , 



(51) 
otherwise , 



with R = 0.3 and r = ^/x^^^y^■ The exact amplitude of the magnetic field 
is arbitrary, since we are solving a linear equation, we used i? = 1. In the 
following, we use exactly the same resolution as Gardiner and Stone (2005). 
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Fig. 5. Magnetic loop advection test: AMR result with the U-MUSCL scheme. The 
two upper plots are for Imax = 7, while the two lower plots are for i^ax = 9. The 
right panels show gray-scale images of the magnetic energy, while the left panels 
show the AMR grid (only "oct" boundaries are shown for clarity, but each oct is in 
fact subdivided into 4 children cells). 

Wc perform the numerical integration of the induction equation up to time 
t = 2 with a Courant factor see (34) is equal to 0.8, for which the magnetic 
loop has evolved twice across the computing box. Our first set of runs use a 
regular Cartesian grid with — 128 and Ny — 64. We test the three different 
schemes, to first order (slope limiters were set to zero) and to second order. 
The aim here is to estimate the numerical diffusion of our various schemes. 

Figure 3 shows gray-scale images of the magnetic energy + By for the 
six runs. Maximum field dissipation occurs at the center and boundaries of 
the loop where the current density is initially singular. Second order schemes 
all give very similar results. At first order, the U-MUSCL scheme performs 
slightly better than the other two, with a more isotropic pattern. To estimate 
more quantitatively the numerical diffusion, we have plotted in Figure 4 the 
total magnetic energy in the computational box as a function of time. Perfect 
advection would have given a constant value of Etot — nR'^. As expected, 
first order schemes arc much more diffusive than the second order ones. All 
the latter give almost identical results, Runge-Kutta being the most diffusive, 
followed by C-MUSCL and then U-MUSCL. At first order, the U-MUSCL 
scheme also appears less diflFusive than the two other schemes. 

We now present the results obtained with our AMR implementation using the 
U-MUSCL scheme (C-MUSCL giving almost identical results). We start with 
a base Cartesian grid with A^a- = 8 and Ny = 4, corresponding to imin = 3. It 
is then adaptively refined up to imax, using the following refinement criterion 
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on the magnetic energy E — + By 

max{\A^E\,\AyE\) 



^ + 0.01 



> 0.05 . (52) 



With this criterion, each cell for which the change of local magnetic energy 
exceeds 5% of the local magnetic energy is refined. The first test is done 
with ijnax = 7, in order to reach the same spatial resolution as the previous 
simulations with a 128 x 64 Cartesian grid. The magnetic energy map att — 2 
is shown in Figure 5, together with a line plot showing the corresponding 
AMR grid. In this last plot, only "oct" boundaries are shown for clarity (each 
oct is in fact subdivided into four children cells). We conclude that the AMR 
results are indistinguishable from the equivalent resolution Cartesian run, but 
the computational cost^ is lower: at time t — 2, the total number of leaf 
cells in the AMR tree is 3149. This is to be compared with the number of 
cells in a Cartesian grid equivalent to the finer resolution which would be 
128 X 64 = 8129. 

In order to illustrate more convincingly the interest of using an AMR grid in 
this case, we have performed the same simulation with now imax = 9. The 
magnetic energy map and the corresponding AMR grid are shown in Figure 5. 
Refinements are now much more localized at the center and boundaries of the 
magnetic loop. Numerical diffusion has dramatically decreased, as shown on 
Figure 4, where the time history of the total magnetic energy is plotted. The 
agreement with the ideal case has improved substantially. The total number of 
cells at t = 2 is now 16433. This is only a factor of 2 greater than the previous 
Cartesian runs, but a factor of 8 lower than the Cartesian grid equivalent to 
the finer resolution 512 x 256 = 131072. 



4-2 The Ponomarenko Dynamo 



One of the simplest known dynamo fiows, and the one we will start our investi- 
gation with, is the Ponomarenko dynamo (Ponomarenko, 1973). The geometry 
of the flow is remarkably simple. In cylindrical polar coordinates (s, 0, z), it is 



V = < 



(0,5^,^^) 





for s < So , 
for s > Sq . 



(53) 



The actual computing time is in our case directly proportionnal to the number 
of active cells. 
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This flow features an abriipt discontinuity across the cyhnder ai s — so, such 
discontinuity yields an intricate behavior in the hmit Rm — > oo. The growth 
rate remains constant in this hmit, but the flow does not quahfy as a proper 
fast dynamo, for the critical eigenmode keeps changing with Rm (see Childress 
and Gilbert, 1995). Variants of this flow, known as "smoothed Ponomarenko 
flows" introduce a typical length scale over which the flow vanishes, and can 
help circumvent this difficulty (Gilbert, 1988). We will however consider here 
the original Ponomarenko flow with an abrupt discontinuity. Since the flow is 
discontinuous, an explicit physical resistivity (associated with a finite value of 
the Reynolds number Rm) is essential in setting the typical lengthscale of the 
magnetic field {£ ~ Rm~^/^). 

As with most dynamo problems, numerical resolution is classically achieved 
using spectral expansions (e.g. Childress and Gilbert, 1995). We use here our 
numerical approach to validate our scheme as well as to test the properties of 
the AMR implementation and its ability to deal with a discontinuous input 
flow. Because of the cylindrical nature of the flow, it is natural to think of 
adapting the scheme to this system of coordinates. We have therefore written 
a cylindrical version of our algorithm (note however that AMR has not been 
implemented in this version of the code). The discontinuity at s — sq corre- 
spond exactly to a cell boundary. It is important to appreciate that there is 
no flow along the s direction with this approach. This implies that numerical 
diffusion vanishes in this direction. It is only nonzero in the and z directions. 
This emphasizes the importance of physical resistivity to obtain meaningful 
results. 

In most practical work, sharp structures in the flow can occur which are not 
necessarily aligned with the grid (see for example the next application). We 
will therefore solve this same dynamo problem using also a Cartesian grid. 
A very large resolution is needed in order to reach a flne discretisation of 
the cyhnder at s = sq (around which the field is localized over a lengthscale 
£ ~ Rm~^/^). This wiU be achieved using our AMR approach. 

The Ponomarenko flow can be investigated analytically (Ponomarenko, 1973). 
Such an analysis reveals that an exponentially growing solution in time can 
be obtained for Rm = Usq/t] > Rmc ~ 17.7 (where U = \JVl'^Sq + u1). This 
is obtained using a spectral expansion of the variables in z and of the form 
exjp{im(f) + ikz). The most unstable mode (at Rm = Rmc) corresponds to 
Uz — l.Sftso, m — 1 and kcSQ — 0.39. For larger magnetic Reynolds number, 
other modes become unstable. 

Using the cylindrical version of our code, we have numerically calculated the 
magnetic energy growth rates for the Ponomarenko flow for a large range of 
magnetic Reynolds numbers, going from Rm — 16.7 to Rm — 2000. 
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Fig. 6. Growth rate for the Ponomarenko dynamo as a function of the magnetic 

Reynolds number. The soUd curve corresponds to the first unstable mode, and the 
dotted line to its harmonic k = 2 x ■ For both modes, the growth rate first 
increases and then decreases with Rm (as expected from analytical linear theory). 
As the Reynolds number increases, a transition occurs from kcio2xkc ■ The star 
symbol Rm = 400 corresponds to the AMR simulation. 

We use So as unit of length, thus senting sq = 1. The grid extends from 0.2 
to 3.5 in radius and the azimuthal coordinate cover the full 27r range. The 
resolution of the grid is {N^., N^, N^) = (64, 50, 64). For the vertical extent of 
the computational domain Lbox, we consider two different cases: case I, for 
which Lbox = 27r/A;g, with kg being 0.39 and case II for which kg = 0.78. Let us 
recall that the classical numerical approach for this problem relies on a Fourier 
expansion in z. In this case, a single mode k is retained in z to enlighten the 
numerical procedure, the optimal value of kc being obtained after optimization. 
Our numerical approach does not allow this sort of mode selection. Instead, 
we can only fix the ^-periodicity of the computational box. In case I, Lbox 
was chosen to match the wavelength of the most unstable mode. However, 
harmonics of the critical mode, being unstable for large Reynolds numbers, 
can also develop in the computational box (as can be seen for example in the 
figure 6.4 of Plunian and Masse, 2002). This is a known issue, which only 
occurs here because the calculation is not restricted to a single mode in z. 

The transition from the first unstable mode to a higher mode in z occurs 
for Reynolds numbers twice critical. We have been able to follow the first 
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Fig. 7. Ponomarenko dynamo with Rm = 400. Left panel: surface of isovalue 
= 10^ for the magnetic energy density at time t = 200. Right panel: mesh 
geometry (for clarity, only "octs" boundaries are displayed here). 

unstable mode to Reynolds number larger than the transition to k = 2 x /Cc by 
carefully selecting the initial condition (and restricting to short enough time 
integrations). We have also turned our attention to the k — 2 x kc instability 
below the transition by studying a computational box of half the standard 
size in the ^-direction. The resulting diagram is presented in figure 6 . 

When Rm = 16.7, the growth rate a of the magnetic energy was found to be 
negative, as expected. For Rm e [17.7,20] a becomes positive in case I and 
the eigenmode corresponds to k — kc ■ When Rm — 20, it is characterized by 
m = 1, k = kc = 0.39 and a = 3.4 x 10~^. This is in very good agreement 
with linear theory (Ponomarenko, 1973). The growth rate obtained for larger 
Rm is represented by the solid line on figure 6. 

In case II, we use a computational domain with half the vertical extend of 
case I. The growing mode has different properties. It is characterized by m = 1 
and k = 2 X kc = 0.78. Its growth rate as a function of Rm is shown on 
figure 6 using the dotted line. The transition between both modes is clear 
near Rm ~ 30. Unless the initial conditions are carefully chosen and the time 
integration is short enough, the mode k — 2xkc will overcome the first critical 
mode for Rm > 30 . 

In order to validate the AMR implementation, we have also performed simula- 
tions on a Cartesian grid with Rm = 400. The size of the box is Lbox = 27r/0.78, 
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similar to case II described above. For this run, we took i„un = 5 and fmax = 8, 
which has yield a maximum of 751360 cells on the grid (this is a factor of 22 
smaller than the number of cells of a 256^ uniform grid). The refinement strat- 
egy was based on the magnitude of the velocity gradient. The growth rate of 
the magnetic energy in this case was measured to be (Jamr = 0.0562 (see the 
star represented in figure 6). This is in very good agreement with the value 
a = 0.0542 obtained with the cylindrical version of the code for the same 
parameter set. 

The structure of the growing eigenmode in this simulation is illustrated in 
figure 7. The left panel represents surfaces of isovalue of the magnetic energy 
density at t = 200 while the structure of the AMR grid is illustrated 

on the right panel. The grid is only refined at the sharp boundary between 
the inner rotating cylinder and the outer motionless medium. This simulation 
demonstrates both the abihty of the scheme to simulate the Ponomarenko 
dynamo using a Cartesian grid and the possibility to handle discontinuities in 
the flow which are not aligned with the grid. 

4.3 The ABC Dynamo 

We now consider another dynamo flow, known as the ABC-flow (for Arnold- 
Beltrami-Childress) . It is defined by a periodic fiow 

u = ^4 (0, sin X, cos x) + B (cos y, 0, sin y) + C (sin z, cos z, 0) . (54) 

We limit our attention here to the classical case of (A : i? : C) = (1 : 1 : 1). 
Let us stress that this test is fully 3D and requires a significant computational 
effort. 

This flow is known as a fast-dynamo: at large, but finite, Rm, eigcnmodcs 
in the form of cigar-shaped structures develop (e.g. Childress and Gilbert, 
1995). They are very localized in space (again i ~ Rm~^^^), therefore con- 
stituting ideal candidates for a investigation using the AMR methodology. 
Traditionally, these problems have been modeled using spectral methods (e.g. 
Galloway and Prisch, 1986). The choice of the velocity profile in the form of 
Fourier modes was largely guided by the underlying numerical method. More 
recently, Archontis et al. (2003) have investigated this fiow using a staggered 
grid and array valued functions. 

We want to emphasize here that because we are now investigating dynamo 
action at large Rm, the stability properties of the Godunov scheme will be 
essential. This will be particularly true using an AMR grid. The refinement 
strategy will ensure that the physical resistivity dominates on the finer grid 
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Fig. 8. Growth rate for the ABC dynamo as a function of the magnetic Reynolds 
number. This diagram agrees remarkably well with the results obtained using a 
spectral description by Galloway and Frisch 1986 (shown as boxes). The star is 
obtained with the AMR implementation. 

which is centered around the cigar shaped magnetic structures (using a thresh- 
old on II B II). Regions relying on a coarser grid, however, will be dominated 
by the numerical resistivity. The properties of the scheme, both in terms of 
stability and of low numerical resistivity are therefore essential ingredients to 
the success of the AMR methodology. 

Dynamo action associated with this flow is not at all trivial. There are at least 
two regions of instability in the parameter space, one for 8.9 < Rm < 17.5 
and a second for Rm > 27 (see Galloway and Frisch, 1986). This second 
instability has been followed up to Rm of a few thousand. We plan to use 
our methodology to investigate higher values of Rm in the near future. This 
intricate behavior of the growth rate with Rm suggests the use of high enough 
values of the magnetic Reynolds number for convergence study. Otherwise, 
an increase of the resisitivity (decrease in Rm) could yield an increase in the 
growth rate by sampling different regions of instability. 

As in the case of the Ponomarenko dynamo, we have calculated the growth 
rate as a function of Rm. The corresponding graph, using a Cartesian grid 
with {Nx, Ny, Nz) = (128, 128, 128) is presented on figure 8 . This diagram is 
in excellent agreement with the spectral results of Galloway and Frisch, 1986, 
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Fig. 9. ABC dynamo investigated with the AMR strategy at Rm = 159. On the 
left panel: surface of isovalue of the magnetic energy density i?^/2 = 3 x 10^^ at 
time t = 80; on the right panel: the AMR mesh geometry (for clarity, only "octs" 
boundaries are displayed here). 

shown in the same figure as squares. 

We now investigate this dynamo using the AMR scheme. We want to stress 
that using AMR without care for such problems is not free of risk, the grid 
being affected by the solution and vice versa. Although for both the advection 
and Ponomarenko tests, the solution has been well captured using straight- 
forward refinement criteria, the situation is more subtle for the ABC flow, for 
which the field generation is not localized. If the strategy is not adequate, 
some regions of the fiow might not be refined as they should be, and thus 
be subject to a large amount of numerical diffusivity. The choice of the opti- 
mal refinement strategy for the ABC flow is beyond the scope of the present 
study. It could for example be based on various flow properties, such as Lia- 
punov exponents, stagnation points, etc, or on various held properties, such 
as gradients, truncation errors, etc. 

As a flrst step, we have used here a criterion based on the magnetic energy 
density which allows the grid to be easily densifled near the cigar-like struc- 
tures: when the local magnetic energy density on level 5, 6, 7... is respectively 
greater than 4, 16, 64... times the mean energy density, new refinements are 
triggered. This strategy is best applied at large Rm for which the magnetic 
structures are well locahzed. We focus here on Rm = 159 (= 1000/27r). 
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Fig. 10. The ABC dynamo is investigated at Rm = 159 with various resolutions. 
The projected magnetic energy density is represented for each run. The convergence 
is demonstrated on the Cartesian grid and the abihty of the AMR grid to capture 
the solution is assessed. 
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The AMR simulation yields a growth rate of 0.052 after 77 hours of wall-time 
computing using 8 processors. It is evolved until t = 80. At that time, the grid 
is composed of 455659 cells. The structure of the eigenmode and the topology 
of the grid are illustrated in figure 9. For comparison, the Cartesian grid 
simulation with 256^ cells yields a growth rate of 0.055 but requires 138 hours 
to evolve the solution only up to i = 46 and using 64 processors! The AMR 
simulation has therefore allowed a gain in memory of a factor of 37, and 
a speed-up of 25 in time. All our computations are compared on figure 10. 
The first four panels show the projected magnetic energy obtained varying 
the resolution from 32^ to 256^. Computations performed with 128^ and 256^ 
cells reveal very little differences and clearly indicate convergence. The two 
bottom snapshots illustrates the structure of the grid in the AMR simulations 
{left panel) and the projected magnetic energy {right panel). There is a good 
agreement between the AMR simulation and the run performed on the 256^ 
grid (about 10%). 



5 Conclusions and perspectives 



We have shown that the Constrained Transport approach for preserving the 
solenoidal character of the magnetic field could be combined with a Godunov 
method, provided a two-dimensional Riemann solver can be used. We have 
further shown how this could be combined with a MUSCL high order scheme. 
We considered three schemes for the predictive step, each with its own merits. 
For a uniform velocity field, these CT schemes are strictly equivalent to well 
known finite volume schemes on the staggered grid. This important result 
provides additional support to the advection properties of the CT framework. 

We have implemented this strategy on a kinematic dynamo problem, for which 
only the induction equation needs to be considered. We have shown that the 
Godunov framework allows an efficient AMR treatment of fast dynamos, by 
ensuring the numerical stability of the scheme in regions solved with a coarse 
grid (for which the effects of the physical diffusion are vanishing) . 

The approach introduced here clearly needs to be adapted to the full set of 
MHD equations, for which solving the Riemann problem is no longer a trivial 
task. This important step raises several additional difficulties and is the object 
of a forthcoming paper (Fromang et al. (2006)). 
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